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A method and means for image processing 

The present invention relates to methods and means for 
5 image processing and, particularly, though not 

exclusively, to the processing of Nuclear Magnetic 
Resonance image data. 

Dynamic imaging, such as dynamic MR imaging, may involve 
10 the acquisition of a time-series of image data sets in 

which each data set is able to represent an image of the 
subject fully in 2 or 3 dimensions at one of a series of 
successive instants (or brief periods) in time. Dynamic 
imaging is often used to record the internal changes in 
15 properties of a stationary sub j ect which are induced by a 
controlled change-inducing influence. An example of this 
is contrast-enhanced medical MR imaging in which a 
"contrast agent " is introduced into a stationary subject 
(e.g. a person) which is detectable as an increase in the 
20 image contrast/brightness of those internal parts of the 
subject in which the agent is located. 

An analysis of the change in image properties (e.g. in 
image contrast/brightness), over time, of chosen fixed 
25 points or static regions within the imaged volume of the 
subject enables an assessment of the properties of the 
chosen fixed points/regions within that imaged volume. 

For example, the image analysis procedure in breast 
30 dynamic imaging typically involves localisation of a 
lesion within an image of the subject containing a 
contrast agent. According to present image analysis 
procedures, the pixel value of each pixel within an image 
of the subject containing a contrast agent (post-contrast 
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images) has subtracted from it the pixel value of the 
corresponding pixel located at the same image point 
within an image of the same subject taken prior to 
introduction of the contrast agent thereto (pre-contrast 
5 image) . The result is a "subtraction image" which 
ideally includes only the effects of the contrast agent 
upon the subject. 



In many situations the image subject (e.g. region of a 
10 body, such as a breast) contains a high content of 
material (e.g. fat) other than the material of interest 
which may obscure relevant image features (e.g. lesion 
shape features) within a post-contrast image. This is 
especially so in breast-scan images . While the technique 
15 of generating subtraction images aims to remove the 
contribution to, or part of, the image arising from fat 
(and other materials to which the use of the contrast 
agent was not directed) , this can never be fully achieved 
due to the fact that materials (especially fat) other 
20 than the material of interest may often also be enhanced 
to a small extent by the presence of a contrast agent. 
Thus, even a subtraction image is likely to suffer the 
obscuring effects discussed above arising from undesired 
image contrast enhancement. 

25 

Analysis of lesion shape within an image and of the 
pattern of enhancement due to contrast agent uptake may 
make a valuable contribution to diagnosis. The analysis 
of lesion shape can only be considered reliable in cases 
30 where the tumour extent is properly segmented within the 
image being studied. The "traditional" approach to this 
problem is manual segmentation of the tumour region 
within subtraction images (e.g. by a radiologist). Even 
if the shape obscuring tissue is not enhancing, it cannot 
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be eliminated using the subtraction technique because of 
noise present in the pre-contrast and/or post-contrast 
image data. As a result the tumour-outlining contour so 
produced is not always well defined. 

5 

Moreover, merely segmenting the outlining image contour 
of a tumour image may not be satisfactory/sufficient if 
the tumour in question has a heterogeneous structure, 
which is not clearly seen on a subtracted image. As a 
10 result the pattern of contrast enhancement is not always 
clearly seen on the subtraction image. Furthermore, a 
subtraction image alone does not always permit one to 
discriminate between tumour and other enhancing tissues, 
e.g. blood vessels. 

15 

The present invention aims to overcome at least some of 
the above deficiencies of the prior art. 

The present invention, at its most general, proposes to 
20 provide a method and/or means enabling one to identify in 
an image, being one of a time-sequence of images 
recording induced changes in pixel values of successive 
images of a subject (such as post-contrast dynamic 
images) , the contribution to the image arising from the 
25 presence within the image subject of a specified 
material/tissue (such non-enhancing and/or low-enhancing 
material, e.g. fat). This identification is proposed to 
be done using a statistical measure derived from the 
dynamic data (e.g. pixel value changes due to contrast 
30 agent uptake) taken from a plurality of the separate 
images forming the time-sequence. 

In the case of a time-sequence of images recording the 
development of contrast-enhancement in a subject, this is 
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preferably done in respect of individual image pixels (or 
"voxels' 7 in a 3D image) from information contained in the- 
pixel intensity data recording the contrast enhancement 
of the subject at a point or pixel location therein. 
5 Thus, the statistical measure is most preferably derived 
according to the variation over time of the pixel 
intensity at an individual image pixel location (i.e. the 
"uptake curves" representing contrast enhancement due to 
contrast agent uptake over time at a given imaged point) . 

10 

The pixel intensity of a given pixel changes over time 
either predominantly in response to a deliberate induced 
effect (e.g. the presence of a contrast enhancing agent) 
at the location within the subject represented by the 

15 pixel, and/or due to image noise. Different 
material/tissue types respond differently to the presence 
of a contrast enhancing agent and the present invention 
proposes a method and means for utilising this fact in 
identifying different materials within the image of a 

20 subject, using a statistical measure of the induced image 
pixel value changes (over time) which has been found to 
be tissue-specific. 

The invention is particularly suited to identifying fat 
25 and/or fatty tissues within the image of a subject and it 
has been found that the statistical measures described 
herein are particularly well suited to identifying the 
present of such tissues specifically. 

30 It is to be understood that the present invention is not 
limited to the processing of Nuclear Magnetic Resonance 
image data and is applicable to the processing of image 
data acquired via other means and in respect of any 
subject- 
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Accordingly, in a first of its aspects, the present 
invention may provide a method of processing a time- 
sequence of separate image data sets which records 
5 induced changes in pixel values of successive images of a 
subject, each set comprising a plurality of image data 
items which each represent the location of an image pixel 
of the image subject according to a common reference 
frame within which the subject is located, the method 
10 including the steps of: 

(a) selecting from each of a plurality of the 
separate image data sets an image data item which 
represents an image pixel located at the same fixed image 
pixel location, thereby to generate a time-domain image 

15 data set preferably containing only image data items 
which represent an image pixel at the same said pixel 
location; 

(b) determining according to a measure (e.g. 
predetermined) of said induced changes or differences as 

20 between the (preferably all of the) pixel values of the 
image data items of the time-domain image data set 
whether the image data items thereof are associated with 
the presence of a specified material (most preferably a 
tissue, e.g. fat) within the image subject. 

25 

There may then follow an additional step (c) , following 
step (b) , of identifying or classifying the image data 
items of the time-domain image data set as being 
unsuitable for use in the generation of an image of the 
30 subject if they are identified in step (b) as being 

associated with the presence of the specified material or 
tissue (e.g. fat) within the image subject and if the 
specified material or tissue is of a type which it is not 
desired to be included within the image of the subject. 
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Otherwise, the image data items in question may be 
identified/classified as being suitable for use in image 
generation of the subject. 



5 Thus, the time-domain image data set contains information 
regarding exclusively the image pixels representing the 
same location within the common reference frame of each 
of the plurality of successive images in the time- 
sequence of images. Most preferably each and all of the 

10 data items within the time-domain image data set are 
employed in determining whether or not the image data 
items thereof correspond with the specified 
material/tissue (e.g. in determining the measure). For 
example, the measure may be a measure of, or be sensitive 

15 to, differences between each pixel value in the time- 
domain image data set as compared to each other such 
pixel value. The result is that the measure is then 
sensitive to the general dispersion of pixel values 
forming the time-domain image data set, and can provide a 

20 measure of how similar/different the image pixel values 
are collectively. 

By analysing the information contained in this dynamic 
time-domain image data set, one is able to derive 

25 information in respect of any one of the plurality of 

separate image data sets within the time-sequence which 
would not be derivable from that one image data set alone 
- namely, whether a pixel value therein results from a 
specified material or tissue (e.g. fat or fatty tissues) 

30 within the image subject, or otherwise. 

The measure of the induced effect may be any one of a 
number of measures, which enable, from an analysis of the 
time-development of the intensity value of a pixel at the 
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same location within each of the separate images in the 
time-sequence, whether that development arose 
predominantly due to a specific material or tissue (e.g. 
fat) . Such measures as would be readily apparent to the 
5 skilled person may be employed according to the present 
invention. However, the measures described herein have 
been found to be specifically suited to identifying the 
presence of fat and fatty tissues as discussed below. 



10 Where the pixel location of an image is non-responsive or 
only weakly responsive to an inducement to changes in its 
pixel value, when the imaged subject is exposed to a 
deliberate change-inducing influence or agent, the value 
of the pixel at that location is likely to be the same 

15 (or similar) in each of the images within the time- 
sequence. Conversely, if the location is strongly 
responsive to such inducements to changes in its pixel 
value, the value of the pixel at that location is likely 
to change significantly over time within the time 

20 sequence. 



Accordingly, the measure is most preferably defined 
according to a measure of the dispersion of, or 
differences between, the values of pixel intensity 

25 associated with the image data items within the time- 
domain image data set. If only a small degree of 
difference (as between values) in such pixel intensity 
values exists when considering all pixel values 
represented within the time-domain image data set, this 

30 amounts to a small dispersion in pixel intensity values 
and suggests the observed dispersion results from the 
pixel values representing tissue/material (e.g. fat) 
within the imaged subject which is non-responsive or only 
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weakly responsive to the inducements to change. The 
image data items of the time-domain image data set may be 
so classified and may be identified as being not suitable 
for use in the generation of an image of the subject. 
5 Conversely, large values of dispersion suggest the 
observed dispersion arises from materials strongly 
responsive to the inducement to change and may be so 
classified. These image data items of the time-domain 
image may also be classified as being suitable for use in 
10 the generation of an image of the subject. 

Preferably, the image data items of the time-domain image 
are determined as being associated with the presence 
within the image subject of the specified material if the 
15 measure exceeds a predetermined threshold value. 

Preferably, the predetermined threshold is a value of the 
measure of dispersion in the values of pixel intensity 
associated with the image data items, within the time- 
domain image data set. 

20 

For example, the measure (e.g. of dispersion) may be 
determined according to the extent to which the pixel 
values represented by the image data items of the time- 
domain image data set differ from equality (i.e. the 

25 extent to which they do not all share the same value) . 
To this extent, the measure (e.g. of dispersion) may be 
determined according to the degree to which a time-domain 
image vector differs from an identity vector in the same 
vector space. Each vector component of the time-domain 

30 image vector is preferably represented by a unique image 
data item of the time-domain image data set. Preferably, 
the array of vector components is arranged in time- 
sequence. All vector components of the identity vector 
share the same value (preferably, the value 1 (one)). 
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The degree of difference (and therefore the measure e.g. 
of dispersion) may be the angle subtended by the time- 
domain image vector with respect to the identity vector. 
Indeed, the measure of dispersion may be determined 
5 according to any suitable property of the time-domain 
vector in vector space. Preferably, the property in 
question is a property of the vector as a whole (e.g. its 
angle relative to another as above) . 



10 For example, step (b.) of the method of the present 
invention in its first aspect may include: 

forming a time-domain image vector wherein each 
image data item of the time-domain image data set 
represents a separate vector component of the time-domain 

15 image vector (e.g. an array of vector components arranged 
in the time-sequence to form the vector) ; and determining 
the measure according to a property of the time-domain 
image vector. 

20 For example, determining the measure could include 

determining the angle (a) subtended by the time-domain 
image vector with respect to the identity vector 
(preferably in the vector space of the time-domain image 
vector); and, employing the angle (a) as the measure. 

25 The predetermined threshold value discussed herein is 
then preferably a predetermined value of the measure. 



Alternatively, the measure e.g. of dispersion in the 
values of pixel intensity associated with the image data 
30 items within the time-domain image data set may be 

determined according to a property or properties of one 
or more of the principal components of the time-domain 
image vector determined according to a principal 
component decomposition of the time-domain image vector. 
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It has been found that the value of the second principal 
component (that is, the principal component associated 
with the second largest eigenvalue of the Principal 
5 Component Analysis (PCA) ) provides a reliable measure of 
the dispersion in the values of pixel intensity 
associated with the image data items within the time- 
domain image data set. Preferably, the predetermined 
threshold is a predetermined value of the second 
10 principal component of the time-domain image vector. 

Preferably, the predetermined threshold value is exceeded 
if the value of the second principal component is greater 
than zero. 

15 For example, step (b) may include: 

forming a time-domain image vector wherein each 
image data item of the time-domain image data set 
represents a separate vector component of (e.g. an array 
of vector components arranged in the time-sequence to 
20 form...) the time-domain image vector; 

determining the second principal component of the 
time-domain image vector; and, 

employing the value of the second principal 
component as the measure. The predetermined threshold 
25 discussed herein value may then preferably be a 
predetermined value of the measure. 

Where step (b) includes the calculation of the angle (a) 
the predetermined threshold value may be deemed to be 
30 exceeded if the angle (a) is less than a threshold angle 
value (a Q ) . The threshold value (a 0 ) may be arbitrarily 
set by the user, but is preferably determined according 
to the distribution of the angular values (a) of the 
time-domain image vectors associated with a plurality of 
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(preferably all of) the image pixel locations represented 
in the time-sequence of data sets. 

The method may include: 
5 repeating step (a) in respect of a plurality of 

image data items thereby to generate a corresponding 
plurality of time-domain image data sets; and 

forming a corresponding plurality of time-domain 
image vectors with each image data item of a given time- 
10 domain image data set representing a separate vector 

component of (e.g. an array of vector components arranged 
in the time-sequence to form...) the given time-domain 
image vector; then 

determining the angle subtended by each of the time- 
15 domain image vectors with respect to the identity vector 
e.g. in the vector space of the time-domain image vector; 
and, 

determining from the distribution of the values of 
the angles (preferably the values of the natural 

20 logarithm of the angles) of all of the time-domain 

vectors the portion of the angular distribution arising 
mainly, predominantly, or substantially only from the 
presence of the specified material within the image 
subject. The predetermined threshold value may then be 

25 deemed to be exceeded if the angle subtended by the time- 
domain image vector falls within, or contributes to, the 
portion of the angular distribution arising from the 
specified tissue/material (e.g. fat). 

30 Preferably, the threshold angle value ( a 0 ) is the angular 
value which demarcates the portion of the angular 
distribution arising mainly, predominantly, or 
substantially only from the specified tissue/material 
from the other portion (s) of the angular distribution. 
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For example, the portion of the angular distribution 
arising mainly, predominantly, or substantially only from 
the specified tissue/material may be determined according 
to a Normal Distribution having distribution parameters 
5 which cause it to most closely correspond with a portion 
of the angular distribution. Demarcation may be 
performed graphically/visually, or may be performed 
analytically and automatically using, for example, the 
well-known ROC technique in a manner such as would be 
10 readily apparent to the skilled person (''Receiver 
Operating Characteristics": see D J Fink and D 
Christiansen, Electronic Engineer's Handbook, McGraw- 
Hill, 3 rd edition (1989), section 25, ppll9-120). 

The method of processing image data, in its first aspect, 
preferably includes the additional step of: 

(d) replacing by a value of zero the pixel value of 
each image data item of each of the plurality of the 
separate image data sets identified as being unsuitable 
for use in the generation of an image of the subject. 
In this way, separate image data sets of the time- 
sequence may be processed to remove tissue/material- 
specific pixels values therefrom leaving unchanged only 
those pixels which have been processed only if they have 
been deemed suitable for use in image formation. 

Furthermore, the removal of noise from the image data may 
be achieved by decomposing according to a principal value 
decomposition the/each time-domain image vector which has 
30 been deemed suitable for use in image formation, and 
discarding from the decomposed vector those principal 
components which are deemed to arise predominantly from 
random causes. 



15 



20 



25 



WO 2005/048161 PCT/GB2004/004760 

13 

For example, the method may also include representing the 
time-domain image vector in terms of a principal 
component decomposition thereof employing all principal 
component vectors and corresponding principal component 
5 values thereof except: those principal component values 
thereof not exceeding a predetermined magnitude; and, 
replacing by a value determined according to the 
principal component decomposition the pixel value of each 
image data item of each of the plurality of the separate 
10 image data sets identified as being suitable for use in 
the generation of an image of the subject. 



However, it has also been found that the first (largest) 
principal component of a time-domain image vector 

15 generally corresponds to the image data common to all 

images of the time-sequence. Thus, in the case where the 
time-sequence records contrast enhancement beginning with 
a pre-contrast image, the first principal component of 
each time-domain image vector of the sequence 

20 approximately represents the information contained within 
the pre-contrast image data set, that being the first set 
within the time-sequence. Subtraction of this 
information from the time-domain image vectors will 
result in the maximisation of the differences between 

25 different time-domain image vectors (e.g. uptake curves) 
by removing the information which is common to all time- 
domain image vectors (e.g. uptake curves). The term 
"subtraction" is to be understand to refer to the 
reconstruction of a vector using its principal components 

30 except (i.e. omitting) its 1 st principal component. 



Thus, the method may also include representing the time- 
domain image vector in terms of a principal component 
decomposition thereof employing all principal component 
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vectors and corresponding principal component values 
thereof except: those principal component values thereof 
not exceeding a predetermined magnitude; and/or the 
largest principal component value thereof; and, 
5 replacing by a value determined according to the 

principal component decomposition the pixel value of each 
image data item of each of the plurality of the separate 
image data sets identified as being suitable for use in 
the generation of an image of the subject. 

10 

Preferably, the method of processing according to the 
invention in its first aspect, includes repeating steps 
(a) and (b) , or steps (a) to (c) , (including any of the 
preferable/alternative features described above) in 
15 respect of pixels at each pixel location represented 
within the time-sequence of separate image data sets. 
Thus, preferably all image pixels are processed as above. 
Of course, the present invention may provide a method of 
forming an image from image data comprised in the time- 
20 sequence of separate image data sets having been 

processed according to the invention in its first aspect 
(including any of the preferable/alternative features 
described above) . 

It is to be understood that the present invention may 
provide means arranged to, or suitable for, the 
implementation of the invention in its first aspect. 
Accordingly, it is to be understood that the present 
invention encompasses such means. 

For example, in a second of its aspects, the present 
invention may provide an image processing means for 
processing a time-sequence of separate image data sets 
which record induced changes in pixel values of 



25 
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successive images of a subject, each set comprising a 
plurality of image data items which each represent the 
location of an image pixel of the image subject according 
to a common reference frame within which the subject is 
5 located,, the image processing means including: 

(a) selection means for selecting from each of a 
plurality of the separate image data sets an image data 
item which represents an image pixel located at the same 
fixed image pixel location, thereby to generate a time- 

10 domain image data set preferably containing only image 
data items which represent an image pixel at the same 
said image pixel location; 

(b) decision means for determining according to a 
measure of said induced changes as between preferably all 

15 of the pixel values of the image data items of the time- 
domain image data set whether the image data items 
thereof are associated with the presence of a specified 
tissue/material within the image subject. 

20 The apparatus may also include: 

(c) identifying means arranged to identify or classify 
the image data items of the time-domain image data set as 
being unsuitable for use in the generation of an image of 
the subject if they are identified by said decision means 

25 (b) as being associated with the presence of the 

specified material within the image subject and if the 
specified material is of a type which it is not desired 
to be included within the image of the subject. The 
identifying means may be arranged to otherwise identify 

30 the image data items in question as being suitable for 
use in image generation of the subject. 

Preferably, the decision means is arranged to determine 
that the image data items of the time-domain image are 
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associated with the specified material if the measure 
exceeds a predetermined threshold value. Preferably, the 
predetermined threshold is a predetermined value of the 
measure of dispersion in the values of pixel intensity 
5 associated with the image data items within the time- 
domain image data set. 



To this extent, the measure e.g. of dispersion may be 
determined according to the degree to which a time-domain 

10 image vector differs from an identity vector in the same 
.vector space, wherein each vector component of the time- 
domain image vector is represented by a unique image data 
item of the time-domain image data set (e.g. the array of 
vector components being arranged in time-sequence) , and 

15 all vector component of the identity vector share the 

same value. The degree of difference (and therefore the 
measure of dispersion) may be the angle subtended by the 
time-domain image vector with respect to the identity 
vector. 

20 

The decision means may include: 

vector means for forming a time-domain image vector 
wherein each image data item of the time-domain image 
data set represents a separate vector component of the 
25 time-domain image vector (e.g. an array of vector 

components arranged in the time-sequence to form the 
vector) ; 

angle determining means for determining the angle 
{a) subtended by the time-domain image vector with 
30 respect to the identity vector e.g. in the vector space 
of the time-domain image vector. 

The decision means is preferably arranged to employ 
the angle (a) as the measure, and the predetermined 
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threshold value is preferably a predetermined value of 
the predetermined measure. 

Alternatively, or additionally, the decision means may 
5 include: 

vector means for forming a time-domain image vector 
wherein each image data item of the time-domain image 
data set represents a separate vector component of (e.g. 
an array of vector components arranged in the time- 
10 sequence) the time-domain image vector; 

principal component means for determining the second 
principal component of the time-domain image vector; 

wherein the decision means is arranged to employ the 
value of the second principal component as the measure. 
15 The predetermined threshold value is then preferably a 
predetermined value of the measure. Preferably, the 
predetermined threshold is exceeded if value of the 
second principal component is greater than zero. 

20 Where the decision means employs the angle {a), the 
predetermined threshold value may be deemed to be 
exceeded if the angle (a) is less than a threshold angle 
value ( a Q ) . 

25 The selection means is preferably arranged to generate in 
respect of a plurality of image data items a 
corresponding plurality of time-domain image data sets; 

the vector means is preferably arranged to form a 
corresponding plurality of time-domain image vectors with 

30 each image data item of a given time-domain image data 

set representing a separate vector component of (e.g. an 
array of vector components arranged in the time-sequence) 
the given time-domain image vector. 
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The angle determining means is preferably arranged to 
determine the angle subtended by each of the time-domain 
image vectors with respect to the identity vector in the 
vector space of the time-domain image vector; and, 
5 the decision means is preferably arranged to 

determine from the distribution of the values of the 
angles (preferably the values of the natural logarithm of 
the angles) of all of the time-domain vectors the portion 
of the angular distribution arising predominantly, 

10 mainly, or substantially only from the presence of the 
specified material within the image subject. The 
predetermined threshold value may then be deemed to be 
exceeded if the angle subtended by the time-domain image 
vector falls within the portion of the angular 

15 distribution arising from the specified material. 

Preferably, the threshold angle value ( a 0 ) is the angular 
value which demarcates the portion of the angular 
distribution arising substantially only from the 
20 specified material from the other portion (s) of the 
angular distribution. 

The decision means may be arranged to determine the 
portion of the angular distribution arising substantially 
25 only from the specified material according to a Normal 
Distribution having distribution parameters which cause 
it to most closely correspond with a portion of the 
angular distribution. 

30 The image processing means of the present invention, in 
its second aspect, preferably includes: 

(d) data modifying means arranged to replace by a 
value of zero the pixel value of each image data item of 
each of the plurality of the separate image data sets 
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identified by the identifying means (c) as being 
unsuitable for use in the generation of an image of the 
subject . 

5 The principal component means is preferably arranged to 
represent the time-domain image vector in terms of a 
principal component decomposition thereof employing all 
principal component vectors and corresponding principal 
component values thereof except: the largest principal 
10 component value thereof; and those principal component 
values thereof not exceeding a predetermined magnitude. 

The image processing means preferably further includes 
data modifying means arranged to replace by a value 
15 determined according to the principal component 

decomposition the pixel value of each image data item of 
each of the plurality of the separate image data sets 
identified as being suitable for use in the generation of 
an image of the subject. 

20 

The image processing means of the present invention, 
according to its second aspect is preferably arranged to 
process each pixel of each of the time-sequence of 
separate image data sets. 

25 

The image processing means preferably includes image 
forming means for forming an image from image data 
comprised in the time-sequence of separate image data 
sets having been processed by the image processing means 
30 according to the invention in its second aspect 
(including any of the aforementioned 
preferable/alternative features thereof) . 
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The invention may provide computer means programmed to 
perform the method according to the invention in its 
first aspect (including any of the aforementioned 
preferable/alternative features thereof) . The invention 
5 may provide a computer program product containing a 

computer program for performing the method according to 
the invention in its first aspect (including any of the 
aforementioned preferable/alternative features thereof) . 
The invention may provide a computer program for 
performing the method according to the invention in its 
first aspect (including any of the aforementioned 
preferable/alternative features thereof) . The invention 
may provide an image generated according to the method of 
the invention in its first aspect (including any of the 
aforementioned preferable/alternative features thereof) . 

The invention shall now be described in terms of the 
following non-limiting examples with reference to the 
accompanying drawings in which: 

Figure 1 illustrates the distribution of the natural 
logarithm of the value of the angle subtended by the 
time-domain image vectors of each image pixel location of 
a time-sequence of images relative to the identity vector 
in the image space of the time-domain image vectors; 

Figure 2 illustrates a decomposition of the 
distribution illustrated in figure 1 into a Normal 
component most closely corresponding with a Normal 
Distribution, and a residual component being the 
difference between the Normal component and the 
distribution of figure 1; 

Figure 3 illustrates a decomposition of the 
distribution illustrated in figure 1 into two fractions 
("fat" and u non fat") discriminated using the value of 
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the second principal component with threshold value equal 
to zero; 

Figure 4 illustrates the effects upon the 
distribution of figure 1 of not removing background noise 
from the distribution (such removal having occurred in 
the distribution of figure 1) ; 

Figure 5(a) illustrates a sequence of post-contrast 
images processed according to the present invention 
(post-contrast time is incrementing from left to right) ; 

Figure 5(b) illustrates a sequence of post-contrast 
subtraction images according to known methods; 

Figure 5(c) illustrates an enlarged first image of 
the sequence illustrated by figure 5(a) 

Figure 6(a) illustrates an image processed according 
to the present invention; 

Figure 6(b) illustrates a post-contrast subtraction 
image according to known methods; 

Figure 7 illustrates a histogram of first principal 
components; 

Figure 8 illustrates a normalised uptake curve, a 
time-development of a first principal component of an 
uptake vector, and a linear uptake model; 

Figure 9a to 9f illustrate a subject image and 
comparison images subject to different image processing; 

Figures 10a to lOd illustrate a subject image and 
comparison images subject to different image processing; 

Figures 11a to 11c illustrate a sequence of 
processed images and a high-resolution comparison image. 

In the following, the invention shall be illustrated in 
terms of the processing of a dynamic image data set 
recording contrast enhancement in a subject containing a 
high fat content. The invention shall be describes in 
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its application as a method and means of identifying fat 
or fatty tissues within an NMR image. 

Consider a time-sequence of M separate image datasets 
5 each of which contains image data elements acquired by 
Magnetic Resonance Imaging of a subject immediately 
before and subsequent to application of an image contrast 
enhancing agent, and representing the imaged subject in 
three dimensions (i.e. a volume image). Each volume 
10 image is represented as a stack of two-dimensional image 
frames with lattice dimensions N x/ N y , N 2 (respectively 
number of columns and rows in a frame, and a stack 
length) , and total number of image pixels (or "voxels" 
for a volume image) is N = N x x N y x N z . 

15 

Let Y, = (y, y iM ), i = 1, N be a time-domain image data 

set (representing a contrast enhancement agent uptake 
vector, i.e. the vector representation of an uptake curve 
of a single voxel) of a single voxel at a fixed image 

20 point common to each of the M separate image data 

sets within the sequence. The time-domain image data set 
is expressed as a vector in M-dimensional Euclidean space 
wherein each image data item of the time-domain image 
data set represents a separate vector component of an 

25 array of vector components arranged in said time-sequence 
to form the time-domain image vector Y, = (>>, y fM ) - 

One may distinguish between non- and low-enhancing 
tissues, such as fat, and the tissues enhancing to a high 
30 extent. 

We assume that each uptake vector belongs to a certain 
class corresponding to a tissue of specific kind. The 



WO 2005/048161 PCT/GB2004/004760 

23 

dispersion/differences in image pixel values of 
corresponded temporal points of the time-domain image 
vectors (uptake vectors) from the same class are only due 
to the random noise. It follows that 

5 

Where A k is a set of uptake vectors that build the k th 
cluster in the Euclidean space within which the time- 
10 domain image vectors are defined. The quantity jii kj is 

the expected value of the j th value (temporal point) of an 
uptake vector that belongs to the k th cluster, and the 
quantity cr is the standard deviation of the random noise 
which may, for simplicity, be assumed to be the same for 
15 all the dynamic images (this assumption is, however, not 
necessary) . 

We shall base our classification on the following model. 
The dispersion of values of the components ( y fJ ) of an 

20 individual time-domain image vector Y f ={y t]9 ... 9 y lM ) (uptake 

vector) can be modelled by two factors: 

(1) contrast enhancement, being induced and due 
to a deliberate effect, increases the 
dispersion of the expected values /J k}9 ... 9 {i kM 

25 of the corresponded cluster; and, 

(2) random noise. 

For a cluster representing low enhancing tissue, the 
30 difference between the temporal points of an uptake curve 
is mostly due to noise and so the dispersion of the 
expectation values Mk t \>—>MkM is sma ll- 
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It follows that for non-enhancing and low-enhancing 
■tissues the dispersion should be less than is observed in 
a. significantly contrast-enhancing tissue, since the 
5 tissue is not subject to a deliberate or significant 
contrast-enhancing effect, or if there is any such effect 
±t relatively small. 

To establish a measure of the enhancement effect, one may 
10 find a measure of the dispersion of fi k{9 ... 9 /2 kM . This may be 

done using the time-domain vector, Y y = (y # y tM ) , and the 

idem (or "identity" ) vector, l IxA/ =(l,...,l) as follows. 

Considering the uptake curves as vectors in M-dimensional 
Euclidean space (as mentioned above) , the smaller the 

15 dispersion of the expectation values Mk t \>~">MkM ' ^ e c l° ser 
to the idem vector will lay the vectors of the cluster 
. Thus, the angle {a) subtended by a given time- 
domain image vector with respect to the idem vector in 
the vector space of the time-domain image vector is 

20 determined and employed as a measure of the degree of 
dispersion present within that vector. This measure 
reflects both the enhancement effect and the random 
noise. The effect of noise should be taken into account. 
One should find a threshold value a 0 . Should an uptake 

25 -vector posses value a > ao, the uptake vector reflects 

significant enhancement. Statistically this is equivalent 
to testing the hypothesis that the vector represents 
tissue, possessing enhancement bigger than the non- or 
low responsive material (e.g. fat) Later we describe how 

30 to determine the threshold analysing the distribution of 
the statistics (a) . 
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Assuming that the values of a are small: 




(1) 



10 



15 



20 



25 



Where y im = -zr^ M y L , is the average of the temporal points of 
the time-domain vector (uptake vector) . 

The angle (a) is determined in respect of all time-domain 
image vectors (uptake vectors) of all image pixel points 
of the image volume thereby to generate a corresponding 
plurality of such angles. It can be straightforward 
proven mathematically that the values of a will be 
approximately log-normally distributed. Each cluster, A k , 
will possess its own distribution, so for the whole 
volume one shall have a mixed distribution. 

One may determine from the distribution of the values of 
the angles of all of said time-domain vectors the portion 
of the angular distribution arising substantially only 
from the presence within the image subject of a specified 
low-enhancing tissue (e.g. fat) . Once this specified 
portion has been determined, one may then determine that 
the contribution to the pixel values of the image data 
items of a given time-domain image vector (uptake vector) 
arising from the specified material exceeds a 
predetermined threshold value if the angle a subtended by 
-that time-domain image vector falls within the specified 
portion of the angular distribution. 

Figure 1 shows an example of the distribution of the 
values of In (a) for a dynamic breast MR measurement (the 
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solid line) . One can see that this is a mixed 
distribution. One can identify a cluster corresponding to 
the small values of a (shown by the fitted Normal 
Distribution (Gaussian) curve - the dashed line) . The 
5 identified cluster is closer to the idem vector in terms 
of a, i.e. possessing small enhancement effect. It also 
represents statistically the most significant and uniform 
fraction of the population. One may conclude that this 
cluster represents tissue not significantly responsive to 

10 contrast agent, such as fat (at least) , within the imaged 
subject. In general, this cluster represents a tissue 
that is not highly responsive to the contrast agent and 
which is not the tissue of interest to which the use of 
the contrast agent was directed. Hereafter, we shall 

15 refer to the identified cluster as a "cluster of fat" to 
simplify the discussion. 



Fitting the normal distribution curve (the dashed line of 
figure 1) to the fat cluster, the whole angular 

20 distribution (dotted curve of figure 2) can be decomposed 
into normal fraction (the solid curve of figure 2) 
representing fat and the residual fraction shown as the 
dashed curve at the right in figure 2. The residual 
fraction is a mixed distribution corresponding to the 

25 enhancing tissues in the imaged subject. The 
significantly enhancing tissues are much more 
heterogeneous than fat, and therefore fall in a number of 
small overlapping clusters. Using the decomposed 
distribution one can find the optimal threshold, <x 0 , 

30 discriminating between normal (fat) and residual 
populations using for example ROC technique [DJ Fink, D. 
Christiansen (editors), Electronics Engineers 1 Handbook, 
McGraw-Hill, 3rd edition, (1989) : section 25 pp. 119- 
120]. 
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Alternatively , one may wish to select a conservative 
threshold that discriminates only fat tissue without risk 
of loss of any valuable information. It is easy to see on 
5 figure 2 that the threshold corresponds to ln(or)«-2.8 . The 
suppression of image effects caused by fat is done by 
replacing by a zero-vector the time-domain image vectors 

(uptake vectors) Y ; = {y it i^^y fM ) in the time-sequence of 
separate image data sets (the matrix D below) which 
10 possess a value of the angle a < oco (i.e. replacement: 
Y,->(0,...,0)) . 



Alternatively, or additionally, the threshold may be 
determined and applied via a decomposition of the time- 
15 domain image vectors via a Principal Component Analysis 
(PCA) as follows. 



Let D be a dynamic image matrix those columns are 
separates time-domain image vectors (uptake vectors) 

20 Y, = y^y fM ) for respective separates image pixel 

locations "i", and each row is formed from the complete 
contents of successive of the M separate datasets of the 
time-sequence, as follows 



25 Z> = 



y\M 



Let Pj=(P\9— 9 Pm) be the j th principal component (PC) 
vector and is the i th column in P. 
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The value 

M 

Zj.i = (Pj)'Yi - j = 1/ .... M ; i = 1, N 

is the value of the j th principle component of i th time- 
domain image vector (uptake vector) . 



The order of the principal component vectors in P depends 
upon the order of the corresponded eigenvalues (X k ) 
thereof, such that X 1 > X 2 > ... > X M - 



10 It has been found that second principal component is 
associated with the statistical measure (a) , and that 
these two quantities possess a strong negative 
correlation. Analysis of regression of the correlation of 
these two quantities has demonstrated that the value of 

15 the second principal component Z 2/ =0 approximates the 

value of the threshold oc 0 . It follows that if Z 2fi > 0 then 
the vector Y± belongs to the "fat cluster"; otherwise the 
vector Yi represents significantly enhancing tissue. (In 
other words, one may say that the above described "fat 

20 cluster" falls on the positive extent of the 2 nd principal 
component.) In effect, the 2 nd principle component is a 
statistical measure of enhancement equivalent to a. It 
should be understood that alternatives to the use of 
Z 2/ =0 as a threshold value (of 2 nd principal component) 

25 could in principle be used. Figure 3 shows distributions 
of In (a) for populations classified as "fat" and "non- 
fat" using the described above PCA approach. One can see 
that the result coincides very well with the 
aforementioned conservative threshold technique 

30 illustrated in figure 2 (the whole distribution and the 

Gaussian distribution are also shown for comparison - not 
normalised) . 



WO 2005/048161 



29 



PCT/GB2004/004760 



It should be noted that if the aforementioned angular 
distribution method (a statistics) is employed, 
background noise (i.e. noise of empty background 
5 surrounding the imaged subject) may obscure the optimal 
threshold, and therefore is preferably removed. 

Figure 4 illustrates the effects upon the distribution of 
figure 1 of not removing background noise (the solid line 

10 of figure 4) from the distribution (such removal having 
occurred in the distribution of figure 1, reproduced as 
the dashed line in figure 4) . However, one can see from 
figure 4 that the background noise overlaps mainly with 
the significantly enhancing tissues. The PCA approach is 

15 robust to the background noise. 

The present invention may also permit as estimation of 
tissue—specific contrast enhancement effects as follows. 

20 One may describe the time-domain image vector (uptake 
vector data) using the following linear model 

Y g =B g b+H g h k + s,, i e A k , £ = 1,2,... (2) 

25 The first term in eq. 2 is a baseline term. The second 
term describes the specific effect of cluster A* , such 
that the vector (B f b + H f h k ) is collinear with the centroid 
vector of the k th cluster, where fi^ = (fi k ,,...,/z A 2 ) is the 
centroid vector whose components are the expectation 

30 values Mk,\> ~ >Mka • The last term, £, , is the residual 
vector. The vectors b,h A are unit vectors. It is also 
required that b±h A J_€ / for all i 9 k . It follows 
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straightforwardly from the definition of the model that 
the residual term contains only random noise. Indeed, the 
random noise is defined as Y,-^. According to the above 
discussion, e, JL ii k and, therefore, contains only the 
5 component of random noise that is orthogonal to \i k . The 
complimentary component of noise (i.e. one collinear with 
]i k ) is spread across the first and second Principal 
Components (largest and next-largest) of the time-domain 
image vector (i.e. the first two terms of equation (3) 
10 below) . 

Existing subtraction imaging techniques define a baseline 
term as Y fX l, where Yi,i is the first, pre-contrast point 

of the uptake vector Yi and 1 is the identity vector. 
15 These existing techniques estimate the enhancement effect 
merely as Y,-7 M 1, thereby leaving the noise "inside" the 
image data. 

Compare the equation (2) with the following PCA 
20 decomposition of Y: 

M 

Y / = H Z njV» ' Pi -L - -L P M (3) 

The PCA decomposition is a generic linear model of the 
25 data where the meaning of each particular term is 
unknown. However, by postulating equivalence of equation 
(2.) and equation (3) one may assume that: 

(1) the first Principal Component of the series in 
equation (3), namely the quantity Z lt ± pi, is a 

30 "baseline" term; 

(2) a few mid-Components describe the specific 
enhancement effect; and, 
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(3) the remaining Principal Components represents the 
noise which may be discarded from the PCA 
representation of Y ± . 

5 Any suitable methods, such as would be readily apparent 
to the skilled person, may be employed to determine the 
most suitable value for the cut-off number m < M being 
the value such that each j th (j > m) principal component 
can be deemed to contain noise and can be "safely" 
10 discarded. 

Thus the specific tissue effect (arising from 
significantly enhancing tissues) can be estimated using 
the following equation 

15 

Y*i = Hi h =2 Z k/i p k (4) 

In effect this is subtraction of the baseline and 
residual terms of equation (2) . Subtraction of this 

20 information from the time-domain image vectors will 
result in minimisation of the redundant information in 
the image vectors such that those vectors, and the 
information within them, therefore, represent 

substantially only the effects specific to tissue 

25 enhancement post-contrast. 

Thus, the method of image processing may also include 
representing the time-domain image vector in terms of a 
principal component decomposition thereof employing all 
30 principal component vectors and corresponding principal 
component values thereof except: those principal 
component values thereof not exceeding a predetermined 
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magnitude; and/or the largest principal component value 
thereof ; and, 

replacing by a value determined according to the 
principal component decomposition of the pixel value of 
5 each image data item of each of the plurality of the 

separate image data sets identified as being suitable for 
use in the generation of an image of the subject. 

The invention, in a preferred embodiment, provides a 
method of processing a time-sequence of separate image 
data sets (each being either 2D or 3D) which records 
development of contrast enhancement in an imaged subject 
into which a contrast-enhancing agent has been 
introduced, the method containing some or all of the 
following steps: 

1. PCA decomposition of the time-domain vector for each 
pixel location within the time-sequence of separate 
image data sets; followed by, 
2 . Suppression of image noise due to the random (and 
low-enhancing) effects (e.g. of fat) within an 
imaged subject as described above; followed by, 
3. Estimation the specific enhancement effect by 
representing each time-domain image vector according 
to equation (4) (and optionally, the data can be de- 
noised removing a few last terms in the PCA series 
after estimation of the noise cut-off term number 
m) . 

It is to be understood that the present invention may 
30 provide means arranged to, or suitable for, the 

implementation of the image processing method of the 
invention. Accordingly, it is to be understood that the 
present invention encompasses such means. 
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The image processing may be performed on any suitable 
image processing means (e.g. a computer work station) 
which preferably includes image forming means for forming 
an image from image data having been processed by the 
5 image processing means according to the invention. 

The invention may provide a computer program product 
(e.g. computer disk or other means for carrying/storing a 
computer program readable by a computer) containing a 
10 computer program for performing the processing method. 
The invention may provide a computer program for 
performing the method according to the invention. 



Figure 5(a) shows the sequence of post-contrast images 
15 resulted from using the suggested approach. The 
background grey colour in each image outlines eliminated 
fatty tissue. The first image of the sequence on the left 
is a first post-injection image . Rest of the images are 
ordered according to their acquisition times from left to 
20 right. The image of the same subject generated according 
to known image subtraction techniques, and produced using 
the same sequence of datasets as that used for figure 
5(a), is shown for comparison on figure 5(b). The time- 
sequence of datasets used for the example consists of 
25 seven image volumes: two pre- and five post-contrast, 
performed after administration of Gd-DTPA. Temporal 
resolution was 90 sec. The data sets were acquired using 
3D Tl-weighted fast spoiled gradient echo sequence, on 
1.5T Siemens scanner. 

30 

One can see granularities in the lesion and differences 
between tissues due to the contrast dynamics illustrated 
in figure 5(a) that are not seen on the corresponded 
subtracted image of figure 5 (b) . Comparison of these 
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characteristics with ones of the anatomically 
identifiable tissues is used for understanding of the 
structure of the lesion. Figure 5(c) shows zoomed first 
image of the sequence of processed post-contrast images 
5 presented ±n figure 5(a). The dark area at the top of 
the lesion (pointed by the white arrow) was later 
identified as a blood vessel. The dark area at the centre 
(black arrow) is a compact region (shape has been 
assessed using an orthogonal display) , and therefore most 
10 likely is part of the tumour that becomes enhanced later 
than the surrounding tissue. 

Figures 6(a) and (b) show another example of an image 
processed according to the present invention (figure 

15 6(a)) and one generated by simple subtraction of datasets 
(figure 6(b)). In figure 6(a) one can see a dark region 
in the middle of the lesion. This non-homogeneity is not 
seen on the corresponding subtracted image of figure 
6(b). The method was initially valuated using 24 studies. 

20 In all the cases the approach has been shown more 
informative that the traditional subtracted image 
technique. Suppression of the ill-effects of fat upon a 
contrast-enhancement image significantly reduces the 
amount of data to be analysed. This results in faster and 

25 more reliable localisation of lesion, and the tumour 
contour is clearly defined. The method provides a 
capability of fast visual and reliable analysis of the 
data that is otherwise available only from analysis of 
individual uptake curves. Visualisation of differences in 

30 the dynamic contrast characteristics of the tissue 

provides a way for discrimination between tumour and 
blood vessels, and identification of granularities. 
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The complete elimination of the fatty tissue surrounding 
the lesion identifies the tumour's contour clearly, and 
provides a straightforward way for segmenting tumour 
extent and extracting features. 

5 

A further example of the present invention shall now be 
described. 

Effective fat suppression in images is of major 
10 importance for dynamic contrast enhanced MRI of the 
breast. Current approaches rely on subtraction, fat 
saturation or water excitation. Selective saturation or 
excitation require high field uniformity that can often 
present a challenge. Subtraction imaging will not 
15 completely remove fat, since fatty tissue is enhancing to 
some small extent. In addition subtraction increases 
noise of resulted difference image, and as a result the 
tumour outline is not always well defined. 

20 The following example of the invention provides a NN fat 
suppression" method based on multivariate statistical 
analysis and modelling of individual voxel (or pixel) 
dynamic contrast agent uptake curves. That is to say, the 
following example method provides a method which enables 

25 one to identify and remove from the image dataset the 
uptake curves that belong to fatty tissue. In addition, 
the method transforms the dynamic data into a new 4D 
dataset (one temporal dimension and three spatial 
dimensions) with maximised contrast between the enhancing 

30 tissues and improved contrast for lesions. The method 
demonstrates capabilities of identification of a lesion' s 
boundary or contour, inhomogeneities, and visualisation 
of spatio-temporal patterns of contrast enhancement on 
Tl-weighited dynamic studies. 
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The method was applied to 16 Gd-DTPA dynamic studies 
acquired on 1.5 Tesla scanners using Tl-weighted fast 
spoiled gradient echo sequence. In all the cases 
5 anatomically correct fat suppression was observed 
comparing the results with corresponded subtracted and 
fat-saturated images . 

The method provides a way of producing dynamic images 
10 combining high spatial and temporal resolution with fat 
suppression. The method is robust to non-uniformity of 
the magnetic field and heterogeneities of tissues. 

Dynamic contrast enhanced MRI (DCE-MRI) is an effective 
15 tool for diagnosis of breast cancer. The diagnosis relies 
on a number of features including among others the 
lesion's shape and spatio-temporal pattern of contrast 
enhancement. Analysis of the shape can only be reliable 
if the lesion's extent is properly segmented. 

20 

Applying specific linear transformation to the uptake 
vectors, the current approach described herein improves 
contrast between tissues with different contrast 
enhancement dynamics. This provides a way of detecting 
25 granularity in a lesion and discriminating between lesion 
and functional tissue. The approach operates on feature 
spaces of unlimited dimensions. 

A dynamic dataset can be considered as a set of uptake 
30 vectors u f = (u iX ,... 9 u iM ) T each of those represents a time 
evolution of signal intensity of a voxel at specific 
spatial location denoted by the index i = l 9 ... 9 N through the 
sequence of M consecutively acquired volumes (i.e. an 
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intensity uptake curve of the / — th voxel) , and ( ) 7 denotes 
transpose operation (vectors will be considered as column 
vectors . ) 

5 We assume that each uptake vector belongs to a certain 
cluster or group of such vectors possessing multivariate 
normal distribution. The members of the same cluster 
possess physiologically the same uptake curve form and 
the observed differences with the cluster are only due to 
10 the random noise. Assuming that values of white noise at 
different time points of the same uptake curve are 
uncorrelated and equivalently distributed, the scatter 
within the clusters is assumed isotropic. Because of 
effects of heterogeneity of the tissues and non- 
15 uniformity of magnetic field, each tissue will be 
represented by a number of such clusters. It is assumed 
that those factors affect the signal in fat by a spatial 
location specific multiplier. Thus the mean vectors of 
the clusters representing the fatty tissue will be 
20 collinear: 

VL k =H kV i Ff k = l,2 3 ... (5) 

There H k is a multiplying factor, and |i F is a normalized 
25 to unit length vector denoting the direction, shared by 
the mean vectors \i k = (ju f} ,... 9 ju iM ) T of the fat clusters 
(hereafter \i F will be referred to as normalized signature 
vector of fat) . If the vector \i F is known (its estimation 
will be discussed below) , the following linear 
30 transformation will create a new feature space where the 
mean vectors of the fat clusters are aligned with the 
zero vector: 
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(6) 



where (ji F ,u,) denotes dot product, and E is the identity 
matrix- Applying expectation operator to eq. 6 and then 
5 substituting \i k according to eq. 5, one can show that the 
expectation of a feature vector v,- is a zero-vector if 
the corresponding uptake vector u ; belongs to one of the 
clusters representing fat. Hereafter the vectors u, will 
be referred to as original uptake vectors, and v,. - as 
10 residual uptake vectors. 

The transform defined by eq. 6 can be interpreted 
geometrically as subtraction from each original uptake 
vector of its projection on the axis defined by the 

15 normalized signature vector of fat ]i r . It follows that 
the transform matrix E-fi^-fi^ is singular and its rank 
(and therefore the dimension of the new feature space) is 
M-l . According to the above discussion, the residual 
uptake vectors of fat will posses a single cluster with 

20 zero mean vector and isotropic scatter. The residual 
uptake vectors of non-fatty tissues will demonstrate 
differences in contrast agent uptake dynamics as compared 
to the fat, and will obviously possess significantly non- 
Gaussian mixture of multivariate distributions. 

25 

To discriminate the feature vectors of fat and non-fatty 
tissues, one may apply principal component analysis 
(PCA) . PCA will search for a new orthogonal basis in the 
space of residual uptake vectors maximizing consecutively 
30 variance of the projections on each of the axes. Since 
the residual uptake vectors of fat have normal 
distribution with zero mean vector and isotropic scatter, 
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their contribution to the total variance of the 
projections will t>e isotropic, and for any axis, the 
projections of the residual vectors of fat will possess 
normal distribution with zero mean. Thus the search will 
5 be driven by variability of the residual uptake vectors 
of the tissues other than fat. The first (i.e. the most 
significant) principal component accounts for the most 
variability in the data, and, therefore, it will be 
maximally associated with the differences between fat and 
10 non-fatty tissues. 

The eigenimage of the first principal component (i.e. an 
image whose voxels/pixels are values of the first 
principal component of the corresponded residual uptake 
15 vectors) will visualise differences in effect of contrast 
enhancement between fat and other tissues, and the fatty 
tissue on this imag-e will be presented as a random noise 
with zero mean. 

20 Fitting the model normal distribution as shown in Fig. 7, 
decomposes the distribution of values of the first 
principal components into a Gaussian fraction modelling 
the cluster of fa.tty tissue, and the residual non- 
Gaussian fraction. Applying receiver operation 

25 characteristic (ROC) analysis, the optimal threshold 
discriminating the fractions is found. Finally, using 
this threshold, a mask of fatty tissue is created and 
used for elimination of fatty tissue from the images by 
performing a masking operation. 

30 

The mask will be a binary image of the same lattice 
dimensions as the original images where i-th voxel will 
be substituted by value 1, if the vector v, was 
classified as representing non-fatty tissue, and zero 
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otherwise. The operation of masking mentioned above is, 
therefore, voxel wise multiplication of an image and the 
mask. The optimal threshold can also be found 
interactively adjusting the threshold and simultaneously 
5 observing that no important details are cut off on the 
eigenimage. 

Fig. 7 illustrates a histogram of the values of first 
principal components of the residual uptake vectors (the 
data is from one of the examples discussed later) . The 
Gauss curve modelling the fat cluster in the feature 
space is fitted varying mean and standard deviation of 
the Gauss curve and the weight of Gaussian fraction in 
the mixture. 

Perfusion of fat is much slower than of the other 
tissues, and typical dynamic measurement shows only the 
beginning of the distribution phase terminating long 
before the equilibrium. At these conditions, we will 
consider enhancement of fat as approximately linear, i.e. 
[i F will be approximated as a linear function of time post 
contrast agent administration: 

ju Fj = i4-(l + r-fj, 7 = 1,...,M (7) 

25 

where r is relative rate of enhancement of signal in fat 
due to the contrast agent, f y - time of acquisition of the 
j-th volume measured since contrast agent administration, 
and A is a constant normalizing the vector to unit 
30 length. The following iterative procedure is used to 
estimate the rate of enhancement r . Applying increasing 
values of r f starting from r=0, the transformation 
defined by eq. 6 is computed followed by PCA of the 
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residual uptake vectors. Then the Gauss model is fitted 
as discussed above. Finally r the accepted value of the 
parameter is one resulting in the Gauss curve with the 
centroid point closest to the zero. 

5 

PCA is known as a good tool to find axes separating 
mixtures of multivariate Gaussian distributions. The PCA 
applied to the uptake vectors before performing the 
transform defined by the eq. 6 will outline the axis ii F . 

10 

It was found that the first (i.e. the most significant) 
principal component of the uptake vectors is very close 
to p F , and can be used as its rough estimate of (Fig. 
8). The most significant of the remaining principal 
15 components (in this case the second) may be used as an 
axis for discriminant analysis. This principal component 
vector is identified as one possessing the smallest angle 
with the identity vector (a vector of equal values) . 

20 The method can be combined with PCA based approached for 
noise reduction by reconstructing the uptake vectors back 
from their principal components values discarding non- 
significant principal components. 

25 Figure 8 illustrates the result of three different 
estimation methods for estimating the normalized 
signature vector of fat: solid line - linear model; 
dotted - first principle component vector of the original 
uptake vectors; dashed - average uptake curve of ROI 

30 selected within fatty tissue (ROI size - 100 voxels), 
this last method is crude and it is shown only for the 
purpose of comparison. 
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The present approach was applied to 16 dynamic studies of 
patients at genetic risk of breast cancer. Each study 
consists of two pre-contrast and at least four post- 
5 contrast datasets acquired using a 3D Tl-weighted fast 
spoiled gradient echo (FLASH) sequence with TE = 10-14 
ms, TR = 4.2-5 ms, and flip angle - 35° on 1 . 5T scanners 
of various models (mainly Siemens Magnetom and GE Signa) . 
The contrast agent was Gd-DTPA administrated by 

10 intravenous bolus injection. Temporal resolution was 90 
sec. in most cases. Image stack plane orientation was 
coronal, with the matrix size 256x128 voxels and 
corresponding spatial resolution 1.33x1.33 mm. Number of 
slices was 64 to 60, and slice thickness 2.5 mm. As part 

15 of the study protocol, 3D Tl-weighted fat-saturated image 
was acquired using similar sequence with higher spatial 
resolution 0.66x0.89 mm (matrix size 512x384). 

The data was previously checked for motion, and the 
20 background was segmented applying threshold on length of 
the uptake vectors. For each case a mask of fat tissue 
was computed as discussed. To verify the correctness of 
the method, the mask was applied to subtraction image, 
and the result was compared with the original (i.e. 
25 unmasked) subtraction image. In all the cases 
anatomically correct fat suppression was observed. The 
normalized signature vector of fat was computed both 
using iterative algorithm and PCA of the original uptake 
vectors. In both cases similar mask images were obtained. 
30 The values of parameter rof eq. 7 were in the range 
0.1-0.15. 

Finally the computed mask was applied to the dataset of 
residual uptake vectors. It was found that the produced 
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images have improved contrast for enhancing lesions 
comparing to the subtracted images. Here we present two 
examples . 

5 Figures 9(a) to 9(f) and figures 10(a) to 10(d) present a 
study of a patient with invasive carcinoma (mixed lobular 
and no specific type (ductal) carcinoma, grade II.) 
Examples using the data of this study were demonstrated 
earlier (Fig. 7 and 8) . 

10 

Fig. 9 (a) illustrates a subtracted image of a dataset 
acquired at 4 50 sec post contrast agent 
administration after applying fat mask and de- 
noising using PCA; Fig. 9(b) illustrates the 

15 original subtracted image; Fig. 9(c) illustrates the 
difference image of (a) and (b) . For further 
comparison there Is illustrated: Fig. 9(d) shows high 
resolution image with fat saturation (spatial 
resolution 0.66x0.89x2.5 mm); Fig9(e) shows an 

20 eigenimage of the 1 st principal component of the 
feature vectors; and Fig9(f) shows a computed binary 
mask of fatty tissue. 

The lesion is shown by arrow on Fig. 9 (a). 

25 

In Fig . 10 (a) - (d) is shown a sequence of post-contrast 
images after transformation according to the eq. 7 and 
fat bit masking Figs . 10 (a, b) , and corresponding sequence 
of subtracted images Figs . 10 (c, d) . Post-administration 
30 time is printed on the images. The lesion is boxed in all 
the images. Comparing the sequences one can see 
anatomically correct fat suppression. The backgrounds of 
the Figs. 10 (a, b) vary since the images contain signed 
values, and the background grey colour corresponds to the 
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value of zero that outlines the suppressed (i.e. masked 
out) fat. One can see granularity of the lesion better 
than it is seen on the corresponded subtracted images. 
Comparison of these characteristics with ones of the 

5 anatomically identifiable tissues is used for 
understanding of the structure of the lesion. The dark 
area at the top of the lesion (pointed by the white 
arrow) was identified as a blood vessel. The dark area at 
the centre (black arrow) is a compact region (shape has 

10 been assessed using an orthogonal display) , and therefore 
most likely is part of trie tumour that becomes enhanced 
later than the periphery, that is seen comparing the 
images 10(a) and 10(b). 

15 Figures 11(a) to 11(c) demonstrate another study of a 
patient with invasive ductal carcinoma (grade II) . 
Comparing to the previous example, this example 
demonstrates denser parenchyma, and therefore of 
interest. Both studies were acquired using Siemens 

20 Magnetom scanner. 

Figs. 11(a) and 11(b) illustrate a sequence of 
processed images. The lesion is in the left breast 
(shown by arrow) . Ring type enhancement is clearly 
25 demonstrated comparing trie images. Figure 11(c) is 
provided for further comparison and was obtained at 
higher resolution Tl weighted image with fat 
saturation. 

30 Elimination of fatty tissue from the dynamic dataset 
significantly reduces the amount of data to be analysed, 
and suggested linear transformation of the dataset 
improves contrast for functional and lesion tissues. This 
may result in faster and more reliable localisation of 
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lesion, and the tumour contour is clearly defined. 
Visualisation of differences in the dynamic contrast 
characteristics of the tissue enables analysis of spatio- 
temporal pattern of contrast enhancement, identification 
5 of lesion 7 s granularity and discrimination between lesion 
and attached healthy tissues (e.g. blood vessels) . In the 
cases when lesion is surrounded by fatty tissue, the 
complete elimination of fat identifies the tumour's 
contour clearly, and provides a straightforward way for 
10 segmenting tumour extent and feature extraction. 

The method provides a classification tool, which is 
robust to tissue heterogeneity and non-uniformity of 
magnetic field. Also the method does not require feature 
15 space dimension reduction (other existing clustering 
approaches operate on 3D feature spaces.) 

The method has been described in terms of Breast MR 
imaging. However, it can be applied to other types of 
20 dynamic MR studies as will be readily apparent to the 
skilled person. Modifications and variants to the above 
embodiments, such as would be readily apparent to the 
skilled person, are encompassed within the scope of the 
present invention. 

25 



